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Abstract.   A  numerical  procedure  is  introduced  to  solve  the 
one-dimensional  equations  of  gas  dynamics  for  a  cylindrically 
or  spherically  symmetric  flow.   The  method  consists  of  a 
judicious  combination  of  Glimm's  method  and  operator  splitting. 
The  method  is  applied  to  the  problem  of  a  converging 
cylindrical  shock. 


Introduction.   The  one-dimensional  equations  for  an  inviscid. 


non-heat  conducting,  radially  symmetric  flow  can  be  written  in 
the  form 


where 


U  = 


m 


Ut  +F(U)^  =  -W(U) 


F(U)  = 


m  /p  +  p 
^m(e+p)/pj 


and  W(U)  =  (a-l^ 


'  m/r     N 

m   /pr 
m(e+p)/pr_, 


(1) 


(2) 


where  p  is  the  density,  u  is  the  velocity,  m  =  pu  is  the  momentum, 
p  is  the  pressure,  e  is  the  energy  per  unit  volume,  t  is  time, 
r  is  the  space  coordinate  of  symmetry,  a  is  a  constant  which  is  2 
for  cylindrical  symmetry  and  3  foi'  spherical  symmetry,  and  the 
subscripts  refer  to  differentiation.   We  may  write 


P   +  1   -2 


pu 


(3) 


where  7  is  the  ratio  of  specific  heats  (a  constant  greater  than  1) . 

There  are  two  major  problems  involved  in  solving  the  system 
(1)  directly.   The  first  is  the  singular  nature  near  the  axis 
(r  =  0),  that  is,  there  are  singular  terms  proportional  to  l/r. 
The  second  problem  is  that  the  momentum  equation  (the  second 
component  equation  of  (l))  cannot  be  put  in  conservation  form. 

These  problems  cause  major  difficulties  near  the  axis. 
These  are  usually  overcome  by  some  ad  hoc  method  such  as 


extrapolation  (Payne  I956).   Another  approach  has  been  to  treat 
this  as  a  problem  in  Cartesian  coordinates  in  two  space  dimensions 
(Lppidus  1971)- 

In  the  method  described  below  both  of  these  problems  have 
been  completely  eliminated.  Thus  there  is  no  need  to  resort  to 
any  trickery  in  order  to  solve  the  system  (l). 

Outline  of  the  Method.   The  first  step  in  the  problem  is  to  use 
the  method  known  as  operator  splitting  to  remove  the  inhomogeneous 
terms  -W(u)  from  the  system  (1).   Thus  v/e  solve  the  system 

Ut+I^Il)r  =  °  (^) 

which  represents  the  one-dimensional  equations  of  gas  dynamics 
in  Cartesian  coordinates. 

The  method  used  to  solve  system  (4)  is  the  random  choice 
m.ethod  introduced  by  Glimm  ( I965 )  and  developed  for  hydro- 
dynamics by  Chorin  (1976).   Details  of  this  method  will  be  given 
in  the  next  section,  for  completeness. 

Once  system  (4)  is  solved,  the  system  of  ordinary  differential 
equations 

Ut  =  -W(U)  (5) 

is  solved,  where  the  solution  of  system  (4)  is  used  to  determine 
the  inhomogeneous  term.  -¥  in  (5)«  There  are  several  reasons  for 
this  approach,  v/hich  viill  be  discussed  in  later  sections. 


Gllimn's  Method.   Consider  the  nonlinear  system  of  equations  (4). 
Divide  time  into  intervals  of  length  At  and  let  Ar  be  the  spatial 
increment.   The  solution  is  to  be  evaluated  at  times  nAt  where  n 
is  a  nonnegative  integer  at  the  spatial  points  iAr,  where 

i  =  0,±1,±2,  ...    and  at  time  (n  +  -^)At  at  (i  +-i)Ar. 

'^n 
The  method  is  a  two  step  method.   Let  u.  approximate 

.     'vn+l/2  ,  ,    1 ,    .    1 ,   .     ,  I  . 

U(iAr,nAt)  and  u.   ,  approximate  U(  (i+ -^jAr,  (n  +  ^jAt )  in  (4). 

'\.n+l/2  .,  . 

To  find  the  solution  u.   ,  .  consider  the  system  (4)  along  v;ith 

the  piecewise  constant  initial  data 


U(r,nAt) 


(i+^)Ar 


^u-^'   .   r  <  (i+i)A: 


(6) 


~n 


This  gives  a  sequence  of  Riemann  problems.   If  At  < 


Ar 


2(|u|+c) 
where  c  is  the  local  sound  speed,  the  waves  generated  by  the 

different  Riemann  problems  v^ill  not  interact.   Hence  the  solution 

v(r, t)  to  the  Riemann  problem  can  be  combined  into  a  single  exact 

solution.   See  Figure  1.   Let  ^   be  an  equidistributed  random 

variable  which  is  given  by  the  Lebesgue  measure  on  the  interval 

[-  -^  ,  ^]  .   Define 

See  Figure  2. 

At  each  time  step,  the  solution  is  approximated  by  a  piece- 
wise  constant  function.   The  solution  is  then  advanced  in  time 
exactly  and  the  new  values  are  sampled.   The  method  depends  on 
solving  the  Riemann  problem  exactly  and  inexpensively. 


Chorin  (1976)  (see  also  Sod  (1976)  and  (1977))  modified  an 
iterative  method  due  to  Godunov  (1959)  which  v/ill  nov;  be  described, 
Consider  the  system  (4)  with  the  initial  data 


U(r,0)  =  ]  {8 


The  solution  at  later  times  looks  like  Figure  3,  where  S-]_  and  Sp 
are  either  a  shock  or  centered  rarefaction  wave.   The  region  S^ 
is  a  steady  state.   The  lines  £-,    and  £^   are  separating  the  states. 
The  slip  line  -^  =   u^  separates  the  region  into  two  parts  with 
possibly  different  values  of  p^,  but  equal  values  of  u^  and  p^. 

Using  this  iterative  method  v/e  first  evaluate  p^  in  the 
state  S^.      Define  the  quantity 

If  the  left  wave  is  a  shock,  using  the  jump  condition  U  [p]  =  [pu], 
we  obtain 

M^  =  P^^^^-U^)  =  p^(u^-Up  ,  (10) 

where  U.  is  the  velocity  of  the  left  shock  and  p^  is  the  density 
in  the  portion  of  S^  adjoining  the  left  shock.   Similarly,  define 
the  quantity 

M.  =  ,/  ,,   .  (11) 


If  the  right  v;ave  is  a  shock,  using  the  jump  condition  U  [p]  =  [pu] 
we  obtain 


^r  =  -Pr^^r-  ^r^  =  -p*(u^-U^)  .  (12) 

where  U  is  the  velocity  of  the  right  shock  and  p^  is  the  density 
in  the  portion  of  S^  adjoining  the  right  shock. 

In  either  case  ((9)  or  (lO)  for  M  and  (ll)  or  (12)  for  M^) 


we  obtain 


where 


M^  =/p7p7  Hp*/p^)  ,  (13) 


(14) 


H^)   =  i  (15) 

J-1. 1. 

2/7  1  -  >^^" 
Upon  eliminaticn  of  u^  from  (9)  and  (11)  we  obtain 

P^   Pr 
P.  =  i   /   "  .  (16) 

Equations  (13),  (l^)^  and  (l6)  represent  three  equations  in  three 
unknoivns  for  which  there  exists  a  real  solution.   Upon  choosing  a 
starting  value  p^  (or  M  and  M  ),  we  iterate  using  these  three 
equations.   Vie  chose  p^  =  -g-^P^+Pr^  *'^°^  details  see  Chorin  (1976) 
or  Sod  (1976)). 

After  p^,  M^  and  M.  have  been  determined  v;e  may  obtain  u^ 
by  eliminating  p^  from  equations  (9)  and  (10), 


^* M^TT •  (1^) 

For  a  discussion  of  the  method  of  choosing  the  random  numbers 
most  efficiently  see  Chorin  (1976). 

Solution  of  the  Ordinary  Differential  Equations.   Once  the  solution 

,,  ,   ~n+l 
of  (4)   u.    is  obtained,  we  have  to  solve  a  system  of  ordinary 

differential  equations  (5).   We  approximate  (5)  by 

n+1   ^^n 

W(u.   )  , 


At       -'-i 
or 

uf^  =  u-  -AtwCuJ^""^)  .  (18) 

This  approximation  (18)  is  the  basic  Cauchy-Euler  scheme  which  is 
just  first  order  accurate.  However,  the  Glimm  scheme  is  at  most 
first  order  accurate  so  there  is  no  reason  for  using  a  high  order 
method  for  solving  the  system  of  ordinary  differential  equations. 
Since  this  system  (5)  is  solved  only  at  interior  points  and 
the  scheme  (18)  does  not  require  values  at  r  =  0,  the  singularity 
at  the  axis  is  eliminated. 

Boundary  Conditions.   Boundary  conditions  need  only  be  applied  to 
the  system  (4)  since  the  system  of  ordinary  differential  equations 
(5)  use  only  interior  points.   So  that  with  the  procedure  described 
by  Chorin  (I976)  the  boundary  condition  at  the  axis  (r  =  O)  is 
readily  handled.   The  boundary  condition  is  imposed  on  the  grid 
point  closest  to  r  =  0,  say  i^^Ar.   A  fake  left  state  is  created 
at  (iQ--g-)Ar  by  setting 


Pig-  1/2   =  %W2 


^0-1/2-  "1^+1/2 


\-^/^  '  W 


In  this  way  the  shock  or  rarefaction  wave  will  reflect  which  on 
the  average  is  exact. 

Application  to  a  Converging  Cylindrical  Shock.   Initially,  a 
cylindrical  diaphragm  of  radius  r„  separates  two  uniform  regions 
of  gas  at  rest  as  in  a  shock  tube  with  the  outer  pressure  and 
density  being  larger  than  the  inner  ones.   After  the  diaphragm 
is  ruptured  (t  >  O),  a  shock  wave  is  created  and  travels  into  the 
low  pressure  region  followed  by  a  contact  discontinuity.   A 
rarefaction  wave  travels  into  the  high  pressure  region.   See 
Figure  4. 

It  is  known  that  a  cylindrical  shock  wave  in  a  compressible 
fluid  increases  in  strength  as  it  converges  toward  the  axis. 
This  c£Ln  be  seen  experimentally  in  Perry  and  Kantrowitz  (1951)' 

In  the  example  given  below  the  pressure  and  the  density  in 
the  inner  region  were  set  equal  to  1.0  and  the  pressure  and  density 
in  the  outer  region  were  set  equal  to  4.0.   This  will  produce  a 
shock  with  initial  strength  of  1.93,  a  contact  discontinuity  and 
a  rarefaction  wave.   We  took  Ar  =  0.01.   The  time  step  At  is 
chosen  so  that  the  Courant-Friedrichs-Lewy  condition  is  satisfied, 
i.e. 


II  <inax  (|u|+c)  , 

where  c  is  the  local  sound  speed. 

In  Figure  5  the  pressure  distribution  is  displayed  at  time 
intervals  of  0.2.   The  shock  appears  as  a  rapid  variation  in  p 
which  is  completely  sharp,  i.e.  the  number  of  zones  over  which 
this  variation  takes  place  is  zero.   As  time  increases  the  shock 
propagates  toward  the  axis.   It  is  observed  that  the  strength  of 
the  shock  increases  with  time.   After  the  passage  of  the  shock, 
the  pressure  behind  the  shock  increases.   When  the  shock  arrives 
at  the  axis  it  is  reflected  and  rises  to  a  large  but  finite  value 
and  a  diverging  shock  appears.   It  is  also  observed  that  the 
pressure  at  a  given  time  behind  the  reflected  shock  decreases 
with  time. 

In  Figure  6  the  velocity  of  the  gas  is  displayed.   The 
behavior  is  similar  to  that  of  the  pressure  except  that  the  con- 
verging shock  decreases  the  velocity  from  zero  to  a  negative 
value.   When  the  shock  is  reflected  from  the  axis,  the  diverging 
has  the  effect  of  producing  a  small  positive  (outward)  velocity. 
As  in  the  case  of  the  pressure  profile,  at  a  given  point  behind 
the  converging  shock  the  velocity  increases  with  time  and  behind 
a  diverging  shock  the  velocity  decreases  with  time. 

The  density  and  energy  profiles  are  displayed  in  Figures  7 
and  8  respectively.   The  basic  properties  of  the  shock  are  similar 
to  those  of  the  pressure  distribution,  except  that  the  rise  in 
density  across  the  shock  is  smaller  due  to  a  temperature  increase. 
In  the  density  and  energy  profiles  a  contact  discontinuity  appears, 


It  is  a  result  of  using  Glimm's  scheme  that  the  contact  dis- 
continuity (as  well  as  the  shock  wave)  is  completely  sharp.   The 
contact  discontinuity  propagates  toward  the  axis  behind  the  con- 
verging shock  cLnd  is  traversed  by  the  reflected  (outgoing)  shock. 

In  Figure  9  the  density  profile  where  the  contact  dis- 
continuity and  the  reflected  shock  wave  have  crossed.   For  a 
polytropic  gas  with  the  same  values  of  7,  higher  sound  speeds 
correspond  to  higher  densities  (Courant  and  Friedrichs,  1948). 
The  interaction  of  a  diverging  shock  wave  and  a  contact  dis- 
continuity propagating  toward  the  axis  results  in  a  reflected 
(converging)  shock  (represented  by  (S)  ),  a  contact  discontinuity 
propaLgating  toward  the  axis  (represented  by  (B)  ),  and  a  transmitted 
(diverging)  shock  (represented  by  (C)  ) . 

Conclusions.   This  method  reduces  the  problem  of  solving  the  one- 
dimensional  equations  of  gas  dynamics  for  a  cylindrically  or 
spherically  symmetric  flow  to  solving  the  one-dimensional  equations 
of  gas  dynamics  in  Cartesian  coordinates  and  a  single  system  of 
ordinary  differential  equations,  by  using  operator  splitting. 

The  equations  of  gas  dynamics  are  solved  using  Glimm's 
method  which  keeps  the  shock  waves  and  contact  discontinuities 
perfectly  sharp.   The  ordinary  differential  equations  are  solved 
using  the  Cauchy-Euler  scheme  at  the  interior  points  only  and  for 
one  time  step.   Thus  the  singular  nature  of  the  original  system 
near  the  axis  is  eliminated.   Since  the  equations  of  gas  dynamics 
are  solved  in  Cartesian  coordinates  the  momentum  equation  can  be 
written  in  conservation  form. 


In  all  our  calculations  there  were  100  spatial  grid  points, 
and  it  tSLkes  about  10.3  seconds  on  a  CDC  760O  to  complete  300 
time  steps. 
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f— r ht=(n+l)k 


t=(n+^2)k 


ht  =  nk 
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Fi^ure    1 


Sequence    of  Riemann  problems    on   grid. 
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((i  +  ^^)h,(n+J5)k) 


I t=(n+i5)k 


t=nk 


{l-h)h 


(l+is)h  (i+l)h  ,(i+/2)h 


Figure   2 
Sampling  procedure    for  Glimm's    scheme, 


15 


Figure  3 
Solution  of  Riemann  problem. 
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Figure  4 
Flow  pattern  for  converging  cylindrical  shock. 


15 


Figure    5 
Pressure    profiles    at    time    intervals    of    0.1. 


16 


0.3 


Figure  6 


Velocity  profiles  at  time  intervals  of  0.1. 
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Figure  7 
Density  profiles  at  time  intervals  of  0.1. 
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Figure 


Energy  profiles  at  time  intervals  of  0.1. 
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Figure  9 

Density  profile  after  interaction  of  diverging 
shock  and  contact  discontinuity. 
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